Figure 5 code
Frances Shepherd
2023-11-14
#if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")
#remotes::install_version("Seurat", version = "4.3.0", lib = "~/bin/seurat_4.3.0")
#install.packages("ggplot2")
#install.packages("dplyr")
#BiocManager::install("glmGamPoi")
#install.packages("sctransform")
#install.packages("purrr")
#install.packages("HGNChelper")
#install.packages("openxlsx")
library(Seurat, lib.loc = "~/bin/seurat_4.3.0")## Attaching SeuratObject## Seurat v4 was just loaded with SeuratObject v5; disabling v5 assays and
## validation routines, and ensuring assays work in strict v3/v4
## compatibility modelibrary(sctransform)
library(ggplot2)
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(HGNChelper)
library(openxlsx)
library(stringr)
library(data.table)##
## Attaching package: 'data.table'## The following objects are masked from 'package:dplyr':
##
## between, first, lastlibrary(purrr)##
## Attaching package: 'purrr'## The following object is masked from 'package:data.table':
##
## transpose#Load control and Cal09 infected gene expression matrices (output by CellRanger)
mock.data <- Read10X(data.dir = "")
cal.data <- Read10X(data.dir = "")mock <- CreateSeuratObject(counts = mock.data,
min.features = 200,
project = "mock")## Warning: Feature names cannot have underscores ('_'), replacing with dashes
## ('-')
## Warning: Feature names cannot have underscores ('_'), replacing with dashes
## ('-')cal <- CreateSeuratObject(counts = cal.data,
min.features = 200,
project = "cal")## Warning: Feature names cannot have underscores ('_'), replacing with dashes
## ('-')
## Warning: Feature names cannot have underscores ('_'), replacing with dashes
## ('-')#Add a column which describes the experimental condition for each dataset
mock[["condition"]] <- "mock"
cal[["condition"]] <- "cal"
#Merge two Seurat objects
merged_seurat <- merge(x = mock,
y = cal,
add.cell.id = c("mock", "cal"))
#Calculate mitochondrial genes
merged_seurat[["percent.mt"]] <- PercentageFeatureSet(merged_seurat, pattern = "^MT-")
#Calculate number of genes detected per UMI (give a sense of data complexity)
merged_seurat[["log10GenesPerUMI"]] <- log10(merged_seurat$nFeature_RNA) / log10(merged_seurat$nCount_RNA)
#Create column that replicates the rownames
merged_seurat[["cells"]] <- rownames(merged_seurat@meta.data)QC
Cells captured per condition:
merged_seurat@meta.data %>%
ggplot(aes(x=condition, fill=condition)) +
geom_bar() +
theme_classic() +
theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)) +
theme(plot.title = element_text(hjust=0.5, face="bold")) +
ggtitle("Number of cells per condition")
UMI counts per cell:
merged_seurat@meta.data %>%
ggplot(aes(color=condition, x=nCount_RNA, fill= condition)) +
geom_density(alpha = 0.2) +
scale_x_log10() +
theme_classic() +
ylab("Cell density") +
geom_vline(xintercept = 500, linetype = "dotted") +
geom_vline(xintercept = 100000, linetype = "dotted") +
ggtitle("Number of UMIs per cell")
Number of genes per cell:
merged_seurat@meta.data %>%
ggplot(aes(color=condition, x=nFeature_RNA, fill= condition)) +
geom_density(alpha = 0.2) +
theme_classic() +
scale_x_log10() +
ggtitle("Genes detected per cell") +
geom_vline(xintercept = 600, linetype = "dotted")
Bimodal distribution is typical of NHBE cells.
Number of genes per UMI:
merged_seurat@meta.data %>%
ggplot(aes(color=condition, x=log10GenesPerUMI, fill= condition)) +
geom_density(alpha = 0.2) +
theme_classic() +
scale_x_log10() +
geom_vline(xintercept = 0.8, linetype = "dotted") +
ggtitle("Average genes per UMI\n(Complexity metric)")
Histogram of the mitochondrial gene counts:
merged_seurat@meta.data %>%
ggplot(aes(color=condition, x=percent.mt, fill=condition)) +
geom_density(alpha = 0.2) +
scale_x_log10() +
theme_classic() +
geom_vline(xintercept = 40, linetype = "dotted") +
ggtitle("Percent mitochondrial gene counts")
Higher mitochondrial gene percentages is also typical of NHBE cells. Assess cutoffs of UMIs between 500 and 100,000 and number of unique features/transcripts at 700:
merged_seurat@meta.data %>%
ggplot(aes(x=nCount_RNA, y=nFeature_RNA, color=percent.mt)) +
geom_point() +
scale_colour_gradient(low = "gray90", high = "black") +
#scale_colour_gradientn(colours = "white", "white", "gray", "black",
# breaks = c(0, ))
stat_smooth(method=lm) +
scale_x_log10() +
scale_y_log10() +
theme_classic() +
geom_vline(xintercept = 500, linetype = "dotted") +
geom_vline(xintercept = 100000, linetype = "dotted") +
geom_hline(yintercept = 600, linetype = "dotted") +
facet_wrap(~condition)## `geom_smooth()` using formula = 'y ~ x'## Warning: The following aesthetics were dropped during statistical transformation: colour
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?
## The following aesthetics were dropped during statistical transformation: colour
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?
After applying above filters, plus a filter for MT% greater than 40%. Look at violin plots of the main QC metrics after filtering:
filtered_merged_seurat <- subset(merged_seurat,
subset = (nCount_RNA > 500) &
(nCount_RNA < 100000) &
(nFeature_RNA > 600) &
(log10GenesPerUMI > 0.8) &
(percent.mt < 40))
VlnPlot(filtered_merged_seurat, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3, group.by = "condition")
Check for cell cycle gene expression
s.genes <- cc.genes$s.genes
g2m.genes <- cc.genes$g2m.genes
seurat_phase <- NormalizeData(filtered_merged_seurat) #Normalize the data
seurat_phase <- FindVariableFeatures(seurat_phase, selection.method = "vst") #Find variable genes
seurat_phase <- ScaleData(seurat_phase) #Scale the data## Centering and scaling data matrixseurat_phase <- RunPCA(seurat_phase) # Perform PCA on the variable genes## PC_ 1
## Positive: S100A2, S100A9, KRT5, KRT17, KRT6A, SERPINB3, LY6D, FABP5, TIMP1, KRT13
## SPINK5, SFN, FRMD6, CCND1, DSG3, AGR2, S100A8, IGFBP6, PHLDB2, TUBB
## KRT15, MMP28, TP63, FLNA, IL1RN, KRT14, TNC, COL17A1, LGALS7, RBP1
## Negative: AGBL4, DNAH11, DCDC1, CFAP54, DNAH3, WDR49, SERPINI2, HYDIN, VWA3A, CFAP299
## CDHR3, DNAH9, DNAH7, LMNTD1, RP1, DNAH12, DNAH6, PACRG, DNAI1, ADGB
## FAM227A, CFAP43, SPEF2, DNAAF1, CFAP61, ARMC3, ZBBX, CFAP46, STK33, TTC29
## PC_ 2
## Positive: KRT5, S100A2, COL17A1, FLNA, IGFBP6, SFN, BCAM, LAMA3, TP63, S100A10
## TENM2, DST, CCDC3, FRMD6, KRT17, C16orf74, CAV1, IGFBP7, ANKH, KRT15
## RBBP8, KRT6A, TUBB6, CALD1, FHOD3, TIMP1, CDH13, MMP28, COL4A6, KRT14
## Negative: CD55, WFDC2, LCN2, BPIFB1, C3, CP, MUC16, PIGR, SLPI, TMC5
## SLC6A14, CXCL17, ABCA13, FAM3D, PLEKHS1, MUC1, ANKRD36C, RIMS1, PPP1R9A, CATSPERB
## FER1L6, FCGBP, MSMB, CAPN13, ARHGEF38, MUC5B, CEACAM6, MUC20, MACC1, S100P
## PC_ 3
## Positive: PARD3B, FTO, RARB, SORBS2, ITGB8, LAMC2, STK38L, FMNL2, GPRC5A, RNF213
## HDAC9, ANTXR1, RIMS2, ITPKC, DST, PTK2B, INPP4B, HMGA2, CD55, TNS3
## CATSPERB, WDPCP, TTC6, KLF7, NEDD9, PLEKHG1, ERBB4, ATP10B, NBEA, B4GALNT3
## Negative: FTH1, DYNLL1, CSTB, IFI27, TUBB4B, S100A9, CAPS, MORN2, DYNLT1, C20orf85
## SCGB1A1, SCGB3A1, SLPI, SERPINB3, LCN2, C9orf24, TUBA1A, CRIP1, PIERCE1, TMEM190
## HSP90AA1, WFDC2, PIFO, TPPP3, SAA1, UFC1, CFAP276, CAPSL, OMG, CETN2
## PC_ 4
## Positive: TOP2A, DIAPH3, MKI67, CENPF, ASPM, RRM2, MELK, ANLN, TPX2, KNL1
## NCAPG, KIF14, KIF18B, GTSE1, PBK, ENSG00000284906, CEP55, CENPE, NCAPH, CDK1
## SPC25, KIF15, CKAP2L, BRCA2, ESCO2, KIF4A, BIRC5, TYMS, NDC80, TACC3
## Negative: IGFBP7, DST, COL17A1, CCND2, TNS3, BCAM, IGFBP3, PTPRQ, PALMD, LOXL4
## RBBP8, SOX6, TENM2, MMP13, ANKH, COL12A1, FYB1, FHOD3, LAMA3, PHLDA1
## COL4A6, PMEPA1, SLC7A2, TINAGL1, ABCC3, CDH13, CLU, ADGRL3, TNC, CCDC3
## PC_ 5
## Positive: SERPINB3, KRT13, KRT4, LY6D, SNX31, SERPINB1, SERPINB13, LYPD3, CLCA4, SPINK5
## PLCE1, PLAUR, RNF168, DSG3, GSTA1, LBH, AQP5, MDM2, ABTB3, GPX2
## CADM1, NR4A3, POU2AF1, SMIM35, SCEL, CCDC60, STRBP, ITPKC, IGSF11, RERG
## Negative: SAA1, IFIT3, SAA2, IFIT1, IFI6, ISG15, IDO1, IFITM1, HLA-DRB1, SLC26A4
## MX2, DDX60L, MX1, MT2A, EPSTI1, IFIT2, HERC6, HLA-DPA1, RSAD2, SAMD9L
## HLA-DPB1, DDX60, USP18, CXCL6, HERC5, SLC5A8, RNF213, HLA-DQA1, HLA-DRA, CXCL5seurat_phase <- CellCycleScoring(seurat_phase,
g2m.features = g2m.genes,
s.features = s.genes)## Warning: The following features are not present in the object: MLF1IP, not
## searching for symbol synonyms## Warning: The following features are not present in the object: FAM64A, HN1, not
## searching for symbol synonymsPCA colored by cell cycle phase:
DimPlot(seurat_phase,
reduction = "pca",
group.by= "Phase")
The left side of the graph is dominated by G1 and G2M cells. I will regress out variation due to cell cycle genes in the next step.
Check for MT gene expression
Do the same PCA for mitochondrial gene expression
#Get distribution of MT gene expression:
summary(seurat_phase@meta.data$percent.mt)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.1665 8.4534 14.0768 14.0511 18.8004 39.9950#Based on these values, create a new variable to bin cells based on level of MT gene expression
seurat_phase@meta.data$mitoFr <- cut(seurat_phase@meta.data$percent.mt,
breaks=c(-Inf, 8.4534, 14.0768, 18.8004, Inf),
labels=c("Low","Medium","Medium high", "High"))Plot PCA and color according to their MT gene expression “bin”.
DimPlot(seurat_phase,
size=10,
reduction = "pca",
group.by = "mitoFr")
There are differences on the PC_1, especially in low MT-expressing cells, so I will regress this variation out as well.
Normalization with scTransform
#Split the object into mock and cal conditions
split_seurat <- SplitObject(seurat_phase, split.by = "condition")
#Perform SCtransformation for the mock dataset first
mock <- SCTransform(split_seurat[["mock"]],
vst.flavor = "v2",
vars.to.regress = c("percent.mt", "S.Score", "G2M.Score"))## vst.flavor='v2' set. Using model with fixed slope and excluding poisson genes.## Calculating cell attributes from input UMI matrix: log_umi## Total Step 1 genes: 16548## Total overdispersed genes: 15496## Excluding 1052 genes from Step 1 because they are not overdispersed.## Variance stabilizing transformation of count matrix of size 17313 by 12610## Model formula is y ~ log_umi## Get Negative Binomial regression parameters per gene## Using 2000 genes, 5000 cells##
|
| | 0%
|
|================== | 25%
|
|=================================== | 50%
|
|==================================================== | 75%
|
|======================================================================| 100%## Setting estimate of 173 genes to inf as theta_mm/theta_mle < 1e-3## # of step1 poisson genes (variance < mean): 0## # of low mean genes (mean < 0.001): 826## Total # of Step1 poisson genes (theta=Inf; variance < mean): 182## Total # of poisson genes (theta=Inf; variance < mean): 1788## Calling offset model for all 1788 poisson genes## Found 197 outliers - those will be ignored in fitting/regularization step## Ignoring theta inf genes## Replacing fit params for 1788 poisson genes by theta=Inf## Setting min_variance based on median UMI: 0.04## Second step: Get residuals using fitted parameters for 17313 genes##
|
| | 0%
|
|== | 3%
|
|==== | 6%
|
|====== | 9%
|
|======== | 11%
|
|========== | 14%
|
|============ | 17%
|
|============== | 20%
|
|================ | 23%
|
|================== | 26%
|
|==================== | 29%
|
|====================== | 31%
|
|======================== | 34%
|
|========================== | 37%
|
|============================ | 40%
|
|============================== | 43%
|
|================================ | 46%
|
|================================== | 49%
|
|==================================== | 51%
|
|====================================== | 54%
|
|======================================== | 57%
|
|========================================== | 60%
|
|============================================ | 63%
|
|============================================== | 66%
|
|================================================ | 69%
|
|================================================== | 71%
|
|==================================================== | 74%
|
|====================================================== | 77%
|
|======================================================== | 80%
|
|========================================================== | 83%
|
|============================================================ | 86%
|
|============================================================== | 89%
|
|================================================================ | 91%
|
|================================================================== | 94%
|
|==================================================================== | 97%
|
|======================================================================| 100%## Computing corrected count matrix for 17313 genes##
|
| | 0%
|
|== | 3%
|
|==== | 6%
|
|====== | 9%
|
|======== | 11%
|
|========== | 14%
|
|============ | 17%
|
|============== | 20%
|
|================ | 23%
|
|================== | 26%
|
|==================== | 29%
|
|====================== | 31%
|
|======================== | 34%
|
|========================== | 37%
|
|============================ | 40%
|
|============================== | 43%
|
|================================ | 46%
|
|================================== | 49%
|
|==================================== | 51%
|
|====================================== | 54%
|
|======================================== | 57%
|
|========================================== | 60%
|
|============================================ | 63%
|
|============================================== | 66%
|
|================================================ | 69%
|
|================================================== | 71%
|
|==================================================== | 74%
|
|====================================================== | 77%
|
|======================================================== | 80%
|
|========================================================== | 83%
|
|============================================================ | 86%
|
|============================================================== | 89%
|
|================================================================ | 91%
|
|================================================================== | 94%
|
|==================================================================== | 97%
|
|======================================================================| 100%## Calculating gene attributes## Wall clock passed: Time difference of 1.691986 mins## Determine variable features## Place corrected count matrix in counts slot## Regressing out percent.mt, S.Score, G2M.Score## Centering data matrix## Set default assay to SCT#Then for the cal-infected dataset
cal <- SCTransform(split_seurat[["cal"]],
vst.flavor = "v2",
vars.to.regress = c("percent.mt", "S.Score", "G2M.Score"))## vst.flavor='v2' set. Using model with fixed slope and excluding poisson genes.## Calculating cell attributes from input UMI matrix: log_umi## Total Step 1 genes: 17139## Total overdispersed genes: 16338## Excluding 801 genes from Step 1 because they are not overdispersed.## Variance stabilizing transformation of count matrix of size 18033 by 17955## Model formula is y ~ log_umi## Get Negative Binomial regression parameters per gene## Using 2000 genes, 5000 cells##
|
| | 0%
|
|================== | 25%
|
|=================================== | 50%
|
|==================================================== | 75%
|
|======================================================================| 100%## Setting estimate of 197 genes to inf as theta_mm/theta_mle < 1e-3## # of step1 poisson genes (variance < mean): 0## # of low mean genes (mean < 0.001): 897## Total # of Step1 poisson genes (theta=Inf; variance < mean): 209## Total # of poisson genes (theta=Inf; variance < mean): 1645## Calling offset model for all 1645 poisson genes## Found 211 outliers - those will be ignored in fitting/regularization step## Ignoring theta inf genes## Replacing fit params for 1645 poisson genes by theta=Inf## Setting min_variance based on median UMI: 0.16## Second step: Get residuals using fitted parameters for 18033 genes##
|
| | 0%
|
|== | 3%
|
|==== | 5%
|
|====== | 8%
|
|======== | 11%
|
|========= | 14%
|
|=========== | 16%
|
|============= | 19%
|
|=============== | 22%
|
|================= | 24%
|
|=================== | 27%
|
|===================== | 30%
|
|======================= | 32%
|
|========================= | 35%
|
|========================== | 38%
|
|============================ | 41%
|
|============================== | 43%
|
|================================ | 46%
|
|================================== | 49%
|
|==================================== | 51%
|
|====================================== | 54%
|
|======================================== | 57%
|
|========================================== | 59%
|
|============================================ | 62%
|
|============================================= | 65%
|
|=============================================== | 68%
|
|================================================= | 70%
|
|=================================================== | 73%
|
|===================================================== | 76%
|
|======================================================= | 78%
|
|========================================================= | 81%
|
|=========================================================== | 84%
|
|============================================================= | 86%
|
|============================================================== | 89%
|
|================================================================ | 92%
|
|================================================================== | 95%
|
|==================================================================== | 97%
|
|======================================================================| 100%## Computing corrected count matrix for 18033 genes##
|
| | 0%
|
|== | 3%
|
|==== | 5%
|
|====== | 8%
|
|======== | 11%
|
|========= | 14%
|
|=========== | 16%
|
|============= | 19%
|
|=============== | 22%
|
|================= | 24%
|
|=================== | 27%
|
|===================== | 30%
|
|======================= | 32%
|
|========================= | 35%
|
|========================== | 38%
|
|============================ | 41%
|
|============================== | 43%
|
|================================ | 46%
|
|================================== | 49%
|
|==================================== | 51%
|
|====================================== | 54%
|
|======================================== | 57%
|
|========================================== | 59%
|
|============================================ | 62%
|
|============================================= | 65%
|
|=============================================== | 68%
|
|================================================= | 70%
|
|=================================================== | 73%
|
|===================================================== | 76%
|
|======================================================= | 78%
|
|========================================================= | 81%
|
|=========================================================== | 84%
|
|============================================================= | 86%
|
|============================================================== | 89%
|
|================================================================ | 92%
|
|================================================================== | 95%
|
|==================================================================== | 97%
|
|======================================================================| 100%## Calculating gene attributes## Wall clock passed: Time difference of 2.328911 mins## Determine variable features## Place corrected count matrix in counts slot## Regressing out percent.mt, S.Score, G2M.Score## Centering data matrix## Set default assay to SCTIntegration
mock_cal <- list(mock = mock, cal = cal)
integ_features <- SelectIntegrationFeatures(object.list = mock_cal,
nfeatures = 3000)
# Prepare the SCT list object for integration
mock_cal <- PrepSCTIntegration(object.list = mock_cal,
anchor.features = integ_features)
#Finally perform the actual integration (uses what's called canonical correlation analysis)
integ_anchors <- FindIntegrationAnchors(object.list = mock_cal,
normalization.method = "SCT",
anchor.features = integ_features)## Finding all pairwise anchors## Running CCA## Merging objects## Finding neighborhoods## Finding anchors## Found 32456 anchors## Filtering anchors## Retained 24496 anchorsseurat_integrated <- IntegrateData(anchorset = integ_anchors,
normalization.method = "SCT")## Merging dataset 1 into 2## Extracting anchors for merged samples## Finding integration vectors## Finding integration vector weights## Integrating data#Perform PCA
seurat_integrated <- RunPCA(object = seurat_integrated)## PC_ 1
## Positive: S100A2, KRT5, S100A9, IGFBP6, KRT17, KRT6A, KRT15, KRT13, SERPINB3, LY6D
## RPS18, RPS12, KRT19, TPT1, RPLP1, KRT14, S100A8, RPS2, FABP5, RPL10
## S100A14, TIMP1, EEF1A1, RPL13, CSTA, SFN, HSPB1, TMSB10, PHLDB2, FGFBP1
## Negative: CFAP299, DNAH12, DNAH7, WDR49, AGBL4, CFAP54, HYDIN, LRRIQ1, DNAH6, DNAH5
## DCDC1, DNAH3, SPAG17, TMEM232, ARMC3, PACRG, DNAH11, CFAP47, DNAH9, NEK10
## CDHR3, RFX3, SPATA17, RP1, CFAP43, SERPINI2, ZBBX, ULK4, STK33, FANK1
## PC_ 2
## Positive: S100A2, IGFBP6, KRT14, KRT15, KRT5, LAMA3, KRT17, IGFBP7, DST, KRT6A
## COL17A1, S100A10, SFN, TIMP1, MT2A, LGALS1, TENM2, CCDC3, CAV1, PHLDB2
## BCAM, TP63, FABP5, RBBP8, FLNA, KRT13, CDH13, FRMD6, NRG1, TINAGL1
## Negative: ANKRD36C, BPIFB1, CD55, PIGR, FCGBP, RIMS1, MUC5B, LCN2, RUNX1, ABCA13
## CP, PLAC8, WFDC2, MACC1, CAPN13, CEACAM6, TMC5, SLC6A14, MUC16, C3
## FER1L6, SLPI, MSMB, PLEKHS1, RARRES1, MUC4, SCGB3A1, GPRC5A, HDAC9, NEBL
## PC_ 3
## Positive: SERPINB3, S100A9, LCN2, LY6D, SLPI, WFDC2, GSTA1, SAA2, SERPINB1, TSPAN1
## CAPS, SCGB3A1, KRT4, SMIM22, MORN2, SERPINB4, AQP5, PI3, CYP4B1, CRIP1
## C9orf24, KRT13, C20orf85, S100A6, CSTA, PIERCE1, UPK1B, PIFO, PSCA, TUBA1A
## Negative: DST, IGFBP6, LAMA3, KRT14, ANKRD36C, COL17A1, KRT15, TENM2, RBBP8, CCDC3
## ITGB1, LGALS1, CAV1, LAMC2, GRIP1, CDH13, TNC, BCAM, TNS3, NRG1
## CBLB, FLNA, KRT17, LAMB3, TP63, ANKH, S100A2, MT2A, TINAGL1, FHOD3
## PC_ 4
## Positive: KRT13, LY6D, KRT4, SERPINB3, SERPINB1, CHST9, DENND2C, CSTA, SPINK5, SERPINB13
## MAPK10, PLAUR, LBH, DSG3, RNF168, AQP3, GSTA1, MYO1E, MTSS1, DSP
## S100A9, CLCA4, PLCE1, MAP3K8, ITPKC, RAB11FIP1, LYPD3, MDM2, ZNF750, SERPINB4
## Negative: SAA2, SAA1, BPIFB1, PIGR, LCN2, FCGBP, RARRES1, SLC26A4, IDO1, MSMB
## CEACAM6, MUC5B, CXCL6, C3, WFDC2, DEFB4A, MT2A, HLA-DRB1, IGFBP6, SLPI
## PROM1, KRT14, IGFBP7, LYPD2, SCGB3A1, HLA-DRA, HLA-DPA1, HLA-DPB1, DST, CXCL5
## PC_ 5
## Positive: OMG, TMEM190, IGFBP7, C9orf24, RRAD, CRIP1, MORN2, GPRC5A, CFAP276, ODF3B
## C20orf85, FAM216B, PIERCE1, MACC1, RIIAD1, PLAC8, MKI67, CETN2, SNTN, WDR38
## NWD1, PIFO, PLAAT2, ANKFN1, PLAUR, TUBA1A, C11orf97, CAPSL, CTXN1, MUC16
## Negative: CDC20B, BPIFB1, CCNO, SCGB3A1, FOXN4, CEP128, SMIM35, ANLN, SCGB1A1, SLPI
## LCN2, PLK4, E2F7, CCDC171, DEUP1, MSMB, CEP112, PLCE1, KIF24, PCM1
## SLC4A8, WFDC2, CROCC2, STRBP, STIL, HES6, SDK1, H2BC11, PIGR, PI3# Run UMAP
seurat_integrated <- RunUMAP(seurat_integrated,
dims = 1:40,
reduction = "pca")## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session## 16:39:26 UMAP embedding parameters a = 0.9922 b = 1.112## Found more than one class "dist" in cache; using the first, from namespace 'spam'## Also defined by 'BiocGenerics'## 16:39:26 Read 30565 rows and found 40 numeric columns## 16:39:26 Using Annoy for neighbor search, n_neighbors = 30## Found more than one class "dist" in cache; using the first, from namespace 'spam'## Also defined by 'BiocGenerics'## 16:39:26 Building Annoy index with metric = cosine, n_trees = 50## 0% 10 20 30 40 50 60 70 80 90 100%## [----|----|----|----|----|----|----|----|----|----|## **************************************************|
## 16:39:30 Writing NN index file to temp file /tmp/Rtmpnxirbj/file2cf0af6f206328
## 16:39:30 Searching Annoy index using 1 thread, search_k = 3000
## 16:39:43 Annoy recall = 100%
## 16:39:45 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 16:39:47 Initializing from normalized Laplacian + noise (using irlba)
## 16:39:48 Commencing optimization for 200 epochs, with 1358062 positive edges
## 16:40:20 Optimization finished# Plot UMAP
DimPlot(seurat_integrated, group.by = "condition") 
Clustering cells
seurat_integrated <- FindNeighbors(object = seurat_integrated,
dims = 1:40)## Computing nearest neighbor graph## Computing SNNseurat_integrated <- FindClusters(object = seurat_integrated,
resolution = 0.2)## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 30565
## Number of edges: 1162079
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9481
## Number of communities: 12
## Elapsed time: 4 seconds#Change factor levels of seurat integrated conditions:
seurat_integrated@meta.data$condition <- factor(seurat_integrated@meta.data$condition, levels = c("mock", "cal"))
DimPlot(seurat_integrated,
reduction = "umap",
label = TRUE,
label.size = 6,
split.by = "condition") &
ggtitle("UMAP of integrated data-by cluster") &
labs(caption = "Resolution = 0.2")
Graphing flu infection
#Reset default assay
DefaultAssay(seurat_integrated) <- "RNA"
#Calculate percent flu genes
seurat_integrated[["percent.flu"]] <- PercentageFeatureSet(seurat_integrated, pattern = "Cal-")
summary(seurat_integrated@meta.data$percent.flu)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00000 0.00000 0.01811 0.10368 0.04241 43.25438#Create bins for low, medium, or highly flu-infected cells.
seurat_integrated@meta.data$flu_fraction <- cut(seurat_integrated@meta.data$percent.flu,
breaks=c(-Inf, 0.07, 1, 10, Inf),
labels=c("Uninfected","Low","Medium", "High"))
mycols <- c("lightgray", "#FFFF99", "orange", "red")
DimPlot(seurat_integrated,
size=10,
reduction = "umap",
group.by = "flu_fraction",
split.by = "condition",
cols = alpha(mycols,0.66)) &
ggtitle("Fraction of flu reads per cell\nCombined analysis") &
labs(caption = "Uninfected: < 0.07% flu\nLow: 0.07-1% flu\nMedium: 1-10% flu\nHigh: >10% flu")
Identifying cell types
Find conserved markers in each cluster:
# Create function to get conserved markers for any given cluster. This works similarly to FindAllMarkers but uses the correct analysis with the integrated dataset.
get_conserved <- function(cluster){
FindConservedMarkers(seurat_integrated,
assay = "RNA",
ident.1 = cluster,
grouping.var = "condition",
only.pos = TRUE) %>%
tibble::rownames_to_column(var = "gene") %>%
cbind(cluster_id = cluster, .)
}
#Get markers for all 12 clusters in the integrated dataset
conserved_markers <- map_dfr(c(0:11), get_conserved)## Testing group mock: (0) vs (3, 1, 2, 6, 4, 7, 5, 10, 11, 8, 9)## Testing group cal: (0) vs (4, 5, 1, 2, 3, 8, 6, 9, 11, 7, 10)## Testing group mock: (1) vs (0, 3, 2, 6, 4, 7, 5, 10, 11, 8, 9)## Testing group cal: (1) vs (4, 5, 0, 2, 3, 8, 6, 9, 11, 7, 10)## Testing group mock: (2) vs (0, 3, 1, 6, 4, 7, 5, 10, 11, 8, 9)## Testing group cal: (2) vs (4, 5, 1, 0, 3, 8, 6, 9, 11, 7, 10)## Testing group mock: (3) vs (0, 1, 2, 6, 4, 7, 5, 10, 11, 8, 9)## Testing group cal: (3) vs (4, 5, 1, 0, 2, 8, 6, 9, 11, 7, 10)## Testing group mock: (4) vs (0, 3, 1, 2, 6, 7, 5, 10, 11, 8, 9)## Testing group cal: (4) vs (5, 1, 0, 2, 3, 8, 6, 9, 11, 7, 10)## Testing group mock: (5) vs (0, 3, 1, 2, 6, 4, 7, 10, 11, 8, 9)## Testing group cal: (5) vs (4, 1, 0, 2, 3, 8, 6, 9, 11, 7, 10)## Testing group mock: (6) vs (0, 3, 1, 2, 4, 7, 5, 10, 11, 8, 9)## Testing group cal: (6) vs (4, 5, 1, 0, 2, 3, 8, 9, 11, 7, 10)## Testing group mock: (7) vs (0, 3, 1, 2, 6, 4, 5, 10, 11, 8, 9)## Testing group cal: (7) vs (4, 5, 1, 0, 2, 3, 8, 6, 9, 11, 10)## Testing group mock: (8) vs (0, 3, 1, 2, 6, 4, 7, 5, 10, 11, 9)## Testing group cal: (8) vs (4, 5, 1, 0, 2, 3, 6, 9, 11, 7, 10)## Testing group mock: (9) vs (0, 3, 1, 2, 6, 4, 7, 5, 10, 11, 8)## Testing group cal: (9) vs (4, 5, 1, 0, 2, 3, 8, 6, 11, 7, 10)## Testing group mock: (10) vs (0, 3, 1, 2, 6, 4, 7, 5, 11, 8, 9)## Testing group cal: (10) vs (4, 5, 1, 0, 2, 3, 8, 6, 9, 11, 7)## Testing group mock: (11) vs (0, 3, 1, 2, 6, 4, 7, 5, 10, 8, 9)## Testing group cal: (11) vs (4, 5, 1, 0, 2, 3, 8, 6, 9, 7, 10)cluster_marker_list <- conserved_markers %>%
filter(minimump_p_val < 0.05) %>%
group_by(cluster_id) %>%
arrange(desc(mock_avg_log2FC), .by_group = TRUE)#Load scripts from SCType github
source("https://raw.githubusercontent.com/IanevskiAleksandr/sc-type/master/R/gene_sets_prepare.R"); source("https://raw.githubusercontent.com/IanevskiAleksandr/sc-type/master/R/sctype_score_.R")
# load auto-detection function
source("https://raw.githubusercontent.com/IanevskiAleksandr/sc-type/master/R/auto_detect_tissue_type.R")
db_ = "https://raw.githubusercontent.com/IanevskiAleksandr/sc-type/master/ScTypeDB_full.xlsx";
#Prepare cell type marker database from scType
tissue = "Lung" # e.g. Immune system,Pancreas,Liver,Eye,Kidney,Brain,Lung,Adrenal,Heart,Intestine,Muscle,Placenta,Spleen,Stomach,Thymus
# prepare gene sets
gs_list = gene_sets_prepare(db_, tissue)
# get cell-type by cell matrix
es.max = sctype_score(scRNAseqData = seurat_integrated[["integrated"]]@scale.data,
scaled = TRUE,
gs = gs_list$gs_positive,
gs2 = gs_list$gs_negative)
# merge by cluster
cL_resutls = do.call("rbind", lapply(unique(seurat_integrated@meta.data$seurat_clusters), function(cl){
es.max.cl = sort(rowSums(es.max[ ,rownames(seurat_integrated@meta.data[seurat_integrated@meta.data$seurat_clusters==cl, ])]), decreasing = !0)
head(data.frame(cluster = cl, type = names(es.max.cl), scores = es.max.cl, ncells = sum(seurat_integrated@meta.data$seurat_clusters==cl)), 10)
}))
sctype_scores = cL_resutls %>% group_by(cluster) %>% top_n(n = 1, wt = scores)
# set low-confident (low ScType score) clusters to "unknown"
sctype_scores$type[as.numeric(as.character(sctype_scores$scores)) < sctype_scores$ncells/4] = "Unknown"
sctype_scores[,1:3]## # A tibble: 12 × 3
## # Groups: cluster [12]
## cluster type scores
## <fct> <chr> <dbl>
## 1 0 Airway epithelial cells 6244.
## 2 3 Airway goblet cells 23243.
## 3 1 Basal cells (Airway progenitor cells) 25648.
## 4 2 Clara cells 5810.
## 5 6 Ciliated cells 28190.
## 6 4 Airway goblet cells 11964.
## 7 7 Ciliated cells 2796.
## 8 5 Ciliated cells 1888.
## 9 10 Ciliated cells 696.
## 10 11 Ionocytes 1286.
## 11 8 Basal cells (Airway progenitor cells) 1535.
## 12 9 Ciliated cells 1557.#Merge cell types with cluster IDs on UMAP
seurat_integrated@meta.data$customclassif = ""
for(j in unique(sctype_scores$cluster)){
cl_type = sctype_scores[sctype_scores$cluster==j,];
seurat_integrated@meta.data$customclassif[seurat_integrated@meta.data$seurat_clusters == j] = as.character(cl_type$type[1])
}
#Rename club cells
seurat_integrated@meta.data <- seurat_integrated@meta.data %>%
dplyr::mutate(across('customclassif', str_replace, 'Clara cells', 'Club cells'),
across('customclassif', str_replace, 'Airway goblet cells', 'Secretory and goblet cells'),
across('customclassif', str_replace, 'Basal cells (Airway progenitor cells)', 'Basal cells'))## Warning: There was 1 warning in `dplyr::mutate()`.
## ℹ In argument: `across("customclassif", str_replace, "Clara cells", "Club
## cells")`.
## Caused by warning:
## ! The `...` argument of `across()` is deprecated as of dplyr 1.1.0.
## Supply arguments directly to `.fns` through an anonymous function instead.
##
## # Previously
## across(a:b, mean, na.rm = TRUE)
##
## # Now
## across(a:b, \(x) mean(x, na.rm = TRUE))Idents(seurat_integrated) <- "customclassif"
#Reorder factor levels
DimPlot(seurat_integrated, reduction = "umap", label = FALSE, repel = TRUE, group.by = 'customclassif', order = rev(c("Ciliated cells", "Secretory and goblet cells", "Basal cells (Airway progenitor cells)", "Airway epithelial cells", "Club cells", "Ionocytes"))) &
ggtitle("Integrated UMAP by cell type") &
scale_color_manual(values = c("#1080FF", "#8000FF", "#FB02FF", "#00C13B", "#48EFDA", "#CC9ECB")) &
NoLegend()
UMAP of the flow markers (except for SiR-Tubulin because there is no specific gene for microtubule production):
FeaturePlot(seurat_integrated, features = c("NGFR", "CEACAM6", "TSPAN8"))
Finding variable features between infected and uninfected cell types
Ciliated cells
Compare infected to uninfected condition:
#Read in ISG list
human.isgs <- read.csv("Human ISG list.csv", header = TRUE)
#Filter ISGs to those present in the seurat dataset
human.isgs$present_in_seurat <- human.isgs$isg_name %in% rownames(seurat_integrated)
present_isgs <- subset(human.isgs, present_in_seurat == "TRUE")
ciliated_flu_markers <- FindMarkers(seurat_integrated, ident.1 = "cal", ident.2 = "mock", group.by = "condition", subset.ident = "Ciliated cells")
#Filter by ISGs
ciliated_flu_isgs <- ciliated_flu_markers %>%
filter(row.names(ciliated_flu_markers) %in% c(human.isgs$isg_name) & p_val_adj < 0.05)Secretory cells
Compare infected vs uninfected:
secretory_flu_markers <- FindMarkers(seurat_integrated, ident.1 = "cal", ident.2 = "mock", group.by = "condition", subset.ident = "Secretory and goblet cells")
#Filter by ISGs
secretory_flu_isgs <- secretory_flu_markers %>%
filter(row.names(secretory_flu_markers) %in% c(human.isgs$isg_name) & p_val_adj < 0.05)Heatmap of ISG expression in mock condition
#Set identities of integrated object to the condition
Idents(seurat_integrated) <- "condition"
#Pull out mock dataset
mock_integrated <- subset(seurat_integrated, idents = c("mock"))
#Scale data in the RNA assays
mock_integrated <- ScaleData(mock_integrated, features = present_isgs$isg_name)## Centering and scaling data matrix#Set identities to cell type
Idents(mock_integrated) <- "customclassif"
#Find markers for each cell type
mock_markers <- FindAllMarkers(mock_integrated)## Calculating cluster Airway epithelial cells## Calculating cluster Secretory and goblet cells## Calculating cluster Basal cells (Airway progenitor cells)## Calculating cluster Club cells## Calculating cluster Ciliated cells## Calculating cluster Ionocytes#Filter markers down to ISGs, and filter to only padj < 0.05
mock_isg_markers_all <- mock_markers %>%
filter(rownames(.) %in% present_isgs$isg_name & p_val_adj < 0.05) %>%
group_by(cluster) %>%
arrange(-avg_log2FC, .by_group = TRUE)
#Top 50 only (use for graphing)
mock_isg_markers <- mock_markers %>%
filter(rownames(.) %in% present_isgs$isg_name & p_val_adj < 0.05) %>%
group_by(cluster) %>%
slice_max(n = 50, order_by = avg_log2FC) %>%
arrange(-avg_log2FC, .by_group = TRUE)
#Reorder the cell type levels so that they will be the same between heat maps
levels(mock_integrated) <- c("Ciliated cells", "Secretory and goblet cells", "Basal cells (Airway progenitor cells)", "Airway epithelial cells", "Club cells", "Ionocytes")
#Heatmap of mock condition, cell type specific ISGs
DoHeatmap(mock_integrated, features = mock_isg_markers$gene, assay = "RNA", group.colors = c("#1080FF", "#8000FF", "#FB02FF", "#00C13B", "#48EFDA", "#CC9ECB")) &
ggtitle("ISG Heatmap of Mock-infected cells\nCombined analysis") &
labs(caption = "Cell type identities were determined with full gene set.\nThe most variable genes per cluster were found with FindAllMarkers.\nCell type specific markers were filtered to ISGs and padj < 0.05.\n50 ISGs per cluster are plotted (or as many as are there if less than 50).")
sessionInfo()## R version 4.3.0 (2023-04-21)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: CentOS Linux 7 (Core)
##
## Matrix products: default
## BLAS/LAPACK: /common/software/install/migrated/openblas/0.3.5_gcc8.2.0_multiarch/lib/libopenblasp-r0.3.5.so; LAPACK version 3.8.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## time zone: America/Chicago
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] purrr_1.0.2 data.table_1.14.8 stringr_1.5.0 openxlsx_4.2.5.2
## [5] HGNChelper_0.8.1 dplyr_1.1.3 ggplot2_3.4.4 sctransform_0.4.1
## [9] SeuratObject_5.0.0 Seurat_4.3.0
##
## loaded via a namespace (and not attached):
## [1] RcppAnnoy_0.0.21 splines_4.3.0
## [3] later_1.3.1 bitops_1.0-7
## [5] tibble_3.2.1 polyclip_1.10-6
## [7] lifecycle_1.0.4 Rdpack_2.6
## [9] globals_0.16.2 lattice_0.21-8
## [11] MASS_7.3-58.4 magrittr_2.0.3
## [13] limma_3.58.1 plotly_4.10.3
## [15] sass_0.4.7 rmarkdown_2.25
## [17] plotrix_3.8-4 jquerylib_0.1.4
## [19] yaml_2.3.7 qqconf_1.3.2
## [21] httpuv_1.6.12 sn_2.1.1
## [23] glmGamPoi_1.14.0 spam_2.10-0
## [25] zip_2.3.0 sp_2.1-1
## [27] spatstat.sparse_3.0-3 reticulate_1.34.0
## [29] cowplot_1.1.1 pbapply_1.7-2
## [31] RColorBrewer_1.1-3 multcomp_1.4-25
## [33] abind_1.4-5 zlibbioc_1.48.0
## [35] Rtsne_0.16 GenomicRanges_1.54.1
## [37] BiocGenerics_0.48.1 RCurl_1.98-1.13
## [39] TH.data_1.1-2 sandwich_3.0-2
## [41] GenomeInfoDbData_1.2.11 IRanges_2.36.0
## [43] S4Vectors_0.40.1 ggrepel_0.9.4
## [45] irlba_2.3.5.1 listenv_0.9.0
## [47] spatstat.utils_3.0-4 TFisher_0.2.0
## [49] goftest_1.2-3 spatstat.random_3.2-1
## [51] fitdistrplus_1.1-11 parallelly_1.36.0
## [53] DelayedMatrixStats_1.24.0 leiden_0.4.3
## [55] codetools_0.2-19 DelayedArray_0.28.0
## [57] tidyselect_1.2.0 farver_2.1.1
## [59] matrixStats_1.1.0 stats4_4.3.0
## [61] spatstat.explore_3.2-5 mathjaxr_1.6-0
## [63] jsonlite_1.8.7 multtest_2.58.0
## [65] ellipsis_0.3.2 progressr_0.14.0
## [67] ggridges_0.5.4 survival_3.5-5
## [69] tools_4.3.0 ica_1.0-3
## [71] Rcpp_1.0.11 glue_1.6.2
## [73] mnormt_2.1.1 gridExtra_2.3
## [75] SparseArray_1.2.2 xfun_0.41
## [77] metap_1.9 mgcv_1.8-42
## [79] MatrixGenerics_1.14.0 GenomeInfoDb_1.38.1
## [81] withr_2.5.2 numDeriv_2016.8-1.1
## [83] fastmap_1.1.1 fansi_1.0.5
## [85] digest_0.6.33 R6_2.5.1
## [87] mime_0.12 colorspace_2.1-0
## [89] scattermore_1.2 tensor_1.5
## [91] spatstat.data_3.0-3 utf8_1.2.4
## [93] tidyr_1.3.0 generics_0.1.3
## [95] httr_1.4.7 htmlwidgets_1.6.2
## [97] S4Arrays_1.2.0 uwot_0.1.16
## [99] pkgconfig_2.0.3 gtable_0.3.4
## [101] lmtest_0.9-40 XVector_0.42.0
## [103] htmltools_0.5.7 dotCall64_1.1-0
## [105] scales_1.2.1 Biobase_2.62.0
## [107] png_0.1-8 knitr_1.45
## [109] rstudioapi_0.15.0 reshape2_1.4.4
## [111] nlme_3.1-162 cachem_1.0.8
## [113] zoo_1.8-12 KernSmooth_2.23-20
## [115] parallel_4.3.0 miniUI_0.1.1.1
## [117] pillar_1.9.0 grid_4.3.0
## [119] vctrs_0.6.4 RANN_2.6.1
## [121] promises_1.2.1 xtable_1.8-4
## [123] cluster_2.1.4 evaluate_0.23
## [125] mvtnorm_1.2-3 cli_3.6.1
## [127] compiler_4.3.0 rlang_1.1.2
## [129] crayon_1.5.2 future.apply_1.11.0
## [131] mutoss_0.1-13 labeling_0.4.3
## [133] plyr_1.8.9 stringi_1.8.1
## [135] viridisLite_0.4.2 deldir_1.0-9
## [137] munsell_0.5.0 lazyeval_0.2.2
## [139] spatstat.geom_3.2-7 Matrix_1.6-2
## [141] patchwork_1.1.3 sparseMatrixStats_1.14.0
## [143] future_1.33.0 statmod_1.5.0
## [145] shiny_1.7.5.1 SummarizedExperiment_1.32.0
## [147] highr_0.10 rbibutils_2.2.16
## [149] ROCR_1.0-11 igraph_1.5.1
## [151] bslib_0.5.1